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Abstract. We examine the applicability and viability of methods to obtain 
knowledge about bound-states from information provided solely in Euclidean 
space. Rudimentary methods can be adequate if one only requires information 
about the ground and first excited state and assumptions made about analytic 
^ ■ properties are valid. However, to obtain information from Schwinger functions 

about higher mass states, something more sophisticated is necessary. A method 
based on the correlator matrix can be dependable when operators are carefully 
^ tuned and errors are small. This method is nevertheless not competitive when 

£h ■ an unambiguous analytic continuation of even a single Schwinger function to 

complex momenta is available. 
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1 Introduction 

The probability measure plays a crucial role in quantum field theory. The simplest 
such measure is the Gaussian distribution 



2t 



which is the fundamental solution of the heat equation and the probability density 
that indicates whether a particle characterised initially by a coordinate q will be 
described by q' after an interval r. This distribution, which is related to the free- 
particle density matrix in statistical mechanics, is a probability measure because 
it is positive definite and normalisable [actually, normalised: lim r ^o £t(q, q') = 
5 3 (q — q')]. In Euclidean space a theory's generating functional can truly be 
expressed in terms of a probability measure and the properties of such measures 
make it likely that the rigorous definition of interacting quantum field theories 
is possible in that case. However, in Minkowski space the probability density 
becomes a probability amplitude, through the appearance of "i" in the exponent, 
and that precludes the formulation of a measure theory. 
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The moments of a probability measure are re-point Schwinger functions. They 
correspond to vacuum expectation values of Euclidean fields and may loosely be 
termed Euclidean space Green functions. When certain conditions are met [1], 
analytic continuation of the Schwinger functions yields the Wightman functions 
and one may prove the reconstruction theorem; namely, the complete content 
of a quantum field theory is recovered from the Wightman functions. 1 This is 
the basis for the contemporary belief that the evaluation of a theory's Schwinger 
functions is equivalent to solving that theory. 

While that may be true in principle one must nonetheless develop practical 
means by which to extract physical information from the Schwinger functions. 
The challenge may be illustrated by noting that complete knowledge of a two- 
point Schwinger function in momentum space corresponds only to direct knowl- 
edge of the expectation value at spacelike momenta. A physical particle pole can 
only appear at timelike momentum. Locating such a pole therefore requires that 
all the conditions be met for a unique analytic continuation. If the Schwinger 
function is known only at a discrete set of points; i.e., on a set of measure zero, 
then that is strictly impossible. In this common instance all that is really possible 
is constrained inference and an estimation of the inherent error. 

For us a context for these observations is the problem of determining the 
spectrum of QCD. Much has been learnt about the two-point Schwinger functions 
for quarks, gluons and ghosts. A consensus has been reached; viz., these functions 
cannot be described by a positive definite spectral density. However, the pointwise 
behaviour of the continuation of these functions to timelike momenta is uncertain 
- see, e.g., Refs. [2, 3, 4, 5, 6, 7] and references therein - and is the subject of 
continuing study, e.g., Ref. [8, 9]. 

These two-point functions appear in the kernels of the Bethe-Salpeter equa- 
tions whose solutions provide information about the properties of colour-singlet 
hadrons. The theory and phenomenology of such applications are reviewed in 
Refs. [10, 11, 12, 13, 14]. The study of ground state mesons in the pseudoscalar 
and vector channels has met with success, as evidenced by recent applications to 
heavy-heavy mesons [15] and vector meson electromagnetic form factors [16]. A 
present challenge is the extension to excited states in these and other channels, 
e.g., Refs. [17, 18, 19, 20], and to channels in which the ground state lies above 

1 GeV, e.g., Refs. [21, 22, 23, 24]. 2 A veracious analytic continuation of two-point 
Schwinger functions becomes increasingly important as the bound-state mass 
becomes larger [29]. The possibility that this continuation might be absent re- 
turns us to the question: can reliable information about bound state properties 
be obtained directly from the Schwinger functions? 

Another context for this question is the numerical simulation of lattice- 
regularised QCD. That approach is grounded on the Euclidean space functional 
integral. Schwinger functions are all that it can directly provide. Hence it can only 
be useful if methods are found so that the question can be answered in the affir- 
mative. It is therefore unsurprising that lattice-QCD practitioners have expended 

1 NB. Minkowski space Green functions are constructed from appropriately time-ordered com- 
binations of Wightman functions. 

2 There is also merit in studying the colour antitriplet quark-quark channel since such diquark 
correlations quite likely play an important role in baryon structure, e.g, [25, 26, 27, 28]. 
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much effort on this problem (see, e.g., Refs. [30, 31] and references therein). Our 
study explores the efficacy of the methods devised for lattice-QCD and whether 
it is worthwhile to adapt them to continuum approaches. 

In Sect. 2 we recapitulate on aspects of Dyson-Schwinger equations (DSEs) 
that are relevant to the study of mesons, focusing in particular on the inhomoge- 
neous Bethe-Salpeter equation for colour-singlet three-point Schwinger functions. 
In Sect. 3 we explore and illustrate the application of a simple momentum-space 
method to the problem of extracting meson masses and vertex residues from 
these Schwinger functions. This is followed by an examination of two straightfor- 
ward configuration-space methods for recovering bound-state information from 
Schwinger functions, and a study of the impact of statistical noise on these proce- 
dures. The picture that emerges is that rudimentary methods are often adequate 
if one only requires information about the ground and first excited state, and 
assumptions made about the Schwinger function's analytic structure are valid. 
However, to obtain information from Schwinger functions about heavier states in 
a given channel, something much more refined is necessary. In Sect. 4 we analyse 
one such method; namely, that based on the correlator matrix [32, 33]. Section 5 
is an epilogue. 



2 Inhomogeneous Bethe-Salpeter Equation 

A calculation of the properties of two-body bound states in quantum field the- 
ory may begin with the Poincare covariant Bethe-Salpeter equation (BSE). The 
inhomogeneous BSE for a pseudoscalar quark-antiquark vertex is 3 



ri(k;P) 



tu 



Z475 y + 



xi(q;P) K%(q,k;P), 



(2) 



where r^(k; P) is the vertex Bethe-Salpeter amplitude, with k the relative and P 
the total momentum; r, . . . ,u represent colour, Dirac and flavour matrix indices; 



xUq;P) = S(q + )ri( q] P)S(q- 



(3) 



q± = q ± .P/2; Z4 is the Lagrangian-mass renormalisation constant, described in 
connection with Eq. (9); and J ? represents a Poincare invariant regularisation of 
the integral, with A the regularisation mass-scale [34, 35]. In Eq. (2), K is the 
fully amputated and renormalised dressed-quark-antiquark scattering kernel and 
S in Eq. (3) is the renormalised dressed-quark propagator. It is notable that the 
product SSK is a renormalisation group invariant. 
The solution of Eq. (2) has the form 



irl{k;P) = - lb [iE h {k-P) 

+ 7 • P F 5 (k;P) + 7 • kk ■ PG 5 (k;P) + k^P u H 5 (k; P)] . 



(4) 



3 We work with two degenerate quark flavours and hence Pauli matrices are suflicient to represent 
the flavour structure. In our Euclidean metric: {7^,7^} = 2<5p„; 7^ = 7 M ; and a-b = X^t=i 
For a timelike vector P p , P 2 < 0. More about the metric conventions can be found in Sect. 2.1 
of Ref. [11]. 
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This is the structure necessary and sufficient to ensure Poincare covariance. It 
follows that quark orbital angular momentum is generally present within pseu- 
doscalar and indeed all bound states. This is quantified, e.g., in Ref. [15]. 

The dressed-quark propagator appearing in the BSE's kernel is determined 
by the renormalised gap equation 

Sip)- 1 = Z 2 (i 7 .p + m bm ) + IJ(p), (5) 

Z(p) = 2k J 9 2 D lll/ (p-q)^- lfl S(q)r^q,p), (6) 

where D^y is the dressed-gluon propagator, r u (q,p) is the dressed-quark-gluon 
vertex, and m bm is the /1-dependent current-quark bare mass. The quark-gluon- 
vertex and quark wave function renormalisation constants, Z\^i^-,A 2 ) respec- 
tively, depend on the renormalisation point, £, the regularisation mass-scale and 
the gauge parameter. 

The gap equation's solution has the form 

s(pyl = z{pXe)^' p+M{p2) ^ (7) 

It is obtained from Eq. (5) augmented by the renormalisation condition 

5 '(p)~ 1 |p2 =? 2 >0 = ij-P + m((), (8) 

where m(C) is the renormalised (running) mass: 

Z 2 (C 2 , A 2 ) m bm (A) = Z 4 (C 2 , A 2 ) m(C) . (9) 

Features of the gap equation and its solutions bear upon the radius of convergence 
for a perturbative expansion in the current-quark mass of physical quantities [36]. 

The solution of the inhomogeneous equation, Eq. (2), exists for all values of 
P 2 , timelike and spacelike, with each bound state exhibited as a pole. This is 
illustrated in Fig. 1, wherein the solution is seen to evolve smoothly with P 2 and 
the pole associated with the pseudoscalar ground state is abundantly clear. 4 

Naturally, the numerical determination of the precise location of the first 
pole (ground state) in E${k 2 = 0;P 2 ) will generally be difficult. The task be- 
comes harder if one seeks in addition to obtain the positions of excited states. 
It is for these reasons that the homogeneous equation [39] is usually used. The 
homogeneous pseudoscalar BSE is obtained from Eq. (2) by omitting the driving 
term; viz., the matrix valued constant (1/2)Z475T J . The equation thus obtained 
defines an eigenvalue problem, with the bound state's mass-squared being the 
eigenvalue and its Bethe-Salpeter amplitude, the eigenvector. As such, the homo- 
geneous equation only has solutions at isolated timelike values of P 2 . However, if 

4 NB. The kernel used for this calculation [37] has hitherto provided only stable bound-states. For 
a resonance the pole is not located on the real axis. Nevertheless, in principle such systems can 
also be handled via the BSE. Herein we will mainly overlook such cases because the associated 
complications are not yet relevant: all contemporary studies employ kernels that omit the 
channels associated with physical decays; and unstable states are still not much studied using 
lattice-QCD. 



M.S. Bhagwat, A. H611, A. Krassnigg, CD. Roberts and S. V. Wright 

200 



160- 



x 
<u 

> 



120- 



80- 



40- 



0.3- 



x 

«0.2H 

> 



0.1H 




— Inverse PS BSE Vertex 

— Inverse Scalar BSE Vertex 



0.2 0.4 0.6 Q.i 
P 2 (GeV 2 ) 



■ - ■ PS BSE Vertex Calculation 
— Scalar BSE Vertex 
Pion Pole 



± 



0.2 



0.4 
P 2 (GeV 2 ) 



0.6 



0.8 



Figure 1. Dimensionless pseudoscalar (PS) and scalar amplitudes, -Es(0;P 2 ) and Es(0; P 2 ), 
respectively, obtained by solving an inhomogeneous BSE of the type in Eq. (2) using the 
renormalisation-group-improved rainbow-ladder truncation introduced in Ref. [37]. (We set 
the model's mass-scale lo = 0.33 GeV and used a current-quark mass m u ,d{C,w) — 3.7 MeV, 
fig = 19 GeV. NB. For lj £ [0.3, 0.5] GeV, ground state observables are constant along a tra- 
jectory ujD = (0.72 GeV) 3 =: m 3 g [38].) The vertical dotted-line indicates the position of the 
ground state 7r mass pole; whereas the vertical solid line marks only the coordinate zero. The 
inset shows 1/-E 5 (0; P 2 ) and 1/E S (0;P 2 ). 



one is employing a framework that can only provide reliable information about 
the form of the Bethe-Salpeter vertex at spacelike momenta; i.e., a framework 
that only provides Schwinger functions, then a scheme must be devised that will 
yield the pole positions from this information alone. 



3 Masses from Schwinger Functions: Simple Methods 

3.1 Momentum space: Fade approximant 
One straightforward approach is to focus on 

^(f 2 ) = Ey i 0;P2) (10 

and locate its zeros. It is plain from Fig. 1 that this method can at least be 
successful for the ground state in the pseudoscalar channel. It is important to 
determine whether the approach is also practical for the determination of some 
properties of excited states. 

For the purpose of exploration and illustration, consider a model for an inho- 
mogeneous vertex (three-point Schwinger function) whose analytic structure is 
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Table 1. Parameters characterising the vertex Ansatz, Eq. (11). They were chosen without prej- 
udice, subject to the constraint in quantum field theory that residues of poles in an observable 
projection of a three-point Schwinger function should alternate in sign [17], and ordered such 
that rrii < rm+i, with i = denoting the ground state and i = 1 the first excited state, etc. 
b = 0.78 is the calculated value of Z^ig, A — 200 GeV) used to obtain the curves in Fig. 1. 



i 


Mass, mi (GeV) 


Residue, en (GeV 2 ) 





0.14 


4.23 


1 


1.06 


-5.6 


2 


1.72 


3.82 


3 


2.05 


-3.45 


4 


2.2 


2.8 



known precisely; namely, 

M-\ 

P 2 + m 2 ' 



i=0 



where: for each i, rrii is the bound state's mass and Oj is the residue of the 
bound state pole in the vertex, which is related to the state's decay constant 
(see Appendix A:); and b is a constant that represents the perturbative back- 
ground that is necessarily dominant at ultraviolet total momenta. The particular 
parameter values we employ are listed in Table 1. This Ansatz provides an infor- 
mation sample that expresses essential qualitative features of true DSE solutions 
for colour-singlet three-point Schwinger functions. 

To proceed we employ a diagonal Pade approximant of order N to analyse 
the information sample generated by Eq. (11); viz., 



, , v C + CiX + . . . + cnx n 

JN{X) = — j ; W 



= P\ (12) 



is used to fit a sample of values of 1/V(P 2 ) defined on a discrete P 2 grid, such 
as would be employed in a numerical solution of the BSE. The known ultraviolet 
behaviour of typical vertices requires a diagonal approximant. 5 There are some 
similarities between this approach and that of the inverse amplitudes method in 
Refs. [40, 41]. 

In a confining theory it is likely that a colour-singlet three-point function 
exhibits a countable infinity of bound state poles. Therefore no finite order ap- 
proximant can be expected to recover all the information contained in that func- 
tion. The vertex model presented in Eq. (11) possesses M bound state poles. It is 
reasonable to expect that an approximant of order N < M can at most provide 
reliable information about the first N — 1 bound states, with the position and 



5 NB. A real-world data sample will exhibit logarithmic evolution beyond the renormalisation 
point. No simple Pade approximant can recover that. However, this is not a problem in practical 
applications because the approximant is never applied on that domain. 



M.S. Bhagwat, A. H611, A. Krassnigg, CD. Roberts and S. V. Wright 



7 



3 




2.5- 



— N = 3 ; Ground State 

— N = 3 ; 1 st excited state 

' - ' N = 3 ; 2 nd excited state 



:<*. * • S»' ■ >. : ^ ■■ *" ■■ — ■ ~ — ■ ' — - 




0.5- 



0, 







0.2 



0.4 



0.6 
(GeV 2 ) 



0.8 



max 



Figure 2. Pole positions (mass values) obtained through a fit of Eq. (12) with N = 3 to data 
for 1/V(P 2 ) generated from Eq. (11) with the parameters listed in Table 1. The coordinate 
fmax is described in the text. Horizontal dotted lines indicate the three lightest masses in Table 
1. The ground state mass (solid line) obtained from the Pade approximant lies exactly on top 
of the dotted line representing the true value. 

residue associated with the iV th pole providing impure information that repre- 
sents a mixture of the remaining M — (N — 1) signals and the continuum. We 
anticipate that this is the pattern of behaviour that will be observed with any 
rank-iV approximation to a Schwinger function. The iV-dependence of the Pade 
fit can provide information on this aspect of the procedure. The domain of space- 
like momenta for which information is available may also affect the reliability 
of bound-state parameters extracted via the fitting procedure. Information on 
this possibility is provided by fitting Eq. (12) to the Ansatz data on a domain 
[0, -Pmaxli an d studying the P^^-dependence of the fit parameters. 

We find that a Pade approximant fitted to 1/V(P 2 ) can accurately recover the 
pole residues and locations associated with the ground and first excited states. 
This is illustrated in Fig. 2, which exhibits the Pj^ ax -dependence of the mass- 
parameters determined via a N = 3 Pade approximant. Plateaux appear for three 
isolated zeros, which is the maximum number possible, and the masses defined by 
these zeros agree very well with the three lightest values in Table 1. This appears 
to suggest that the procedure has performed better than anticipated. However, 
that inference is seen to be false in Fig. 3, which depicts the P^^-dependence 
of the pole residues. While the results for oo,i are correct, the result inferred 
from the plateau for is incorrect. Plainly, if the value of 02 were not known 
a priori, then one would likely have been misled by the appearance of a plateau 
and produced an erroneous prediction from the fit to numerical data. A N = 3 
approximant can truly at most only provide reliable information for the first 
N — 1 = 2 bound states. 
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Figure 3. Pole residues obtained through a fit of Eq. (12) with N = 3 to data for l/V(P 2 ) 
generated from Eq. (11) with the parameters in Table 1. Horizontal dotted lines indicate the 
residues associated with the three lightest masses in Table 1. The residue associated with the 
ground state (solid line) lies exactly atop the dotted line representing the true value. The 
residue for the second pole exhibits a plateau at the correct (negative) value. However, the 
plateau exhibited by the result for the residue of the third pole is wrong. 

3.2 A realistic test 

Following this exploration of the viability in principle of using spacelike data 
and the Pade method to extract bound state information, we considered the 
DSE-calculated pseudoscalar vertex that is in depicted Fig. 1. The masses and 
residues (see the appendix) for the ground state pseudoscalar and its first radial 
excitation were obtained from the homogeneous BSE in Ref. [17]. The comparison 
between these masses and those inferred from the Pade approximant is presented 
in Fig. 4. It appears that with perfect (effectively noiseless) spacelike data at 
hand, reliable information on the masses can be obtained. However, while the 
ground state residue is correctly determined, that of the first excited state is not. 
We will subsequently explain this result. 

In the rainbow-ladder truncation of QCD's DSEs the next heaviest ground- 
state meson, after the pseudoscalar, is the scalar [11]. The inhomogeneous scalar 
BSE vertex can be calculated [42] and the analysis described above repeated. The 
relative magnitude of the scalar vertex compared with the pseudoscalar is shown 
in Fig. 1. In the vicinity of P 2 = the signal for a bound state is much suppressed 
because in this model the mass-squared of the quark-core in the ground-state 
scalar meson is P 2 = -(0.675 GeV) 2 = -0.456 GeV 2 . The contribution to the 
signal from a more massive state will be further suppressed. Figure 5 shows the 
results from an N = 3 fit of Eq. (12) to 1/Es(k 2 = 0;P 2 ) calculated using the 
inhomogeneous scalar BSE. The two plateaux give masses in agreement with 
results obtained directly from the homogeneous BSE; i.e., m 1 3i = 0.675 GeV 
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Figure 4. Mass and residue obtained from a fit of Eq. (12) with N = 3 to the DSE result for 
XjE^ik 2 = 0; P 2 ) in Fig. 1. Horizontal dotted lines indicate the masses and residues obtained 
for the ground and first excited state via a direct solution of the homogeneous BSE. 
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Figure 5. Mass and residue obtained from a fit of Eq. (12) with N — 3 to the DSE result for 
l/Es(k 2 = 0;P 2 ); i.e., for the scalar channel, also depicted in Fig. 1. As in Fig. 4, horizontal 
dotted lines indicate the masses and residues obtained for the ground and first excited state via 
a direct solution of the homogeneous BSE. 



and m 2 3L = 1-16 GeV. However, as with the pseudo-scalar, the residue associated 
with the first excited state is wrong. 

The pattern exposed here is not in accord with the expectations raised in 
Sect. 3.1. However, this is only evident because we can simultaneously use time- 
like information to calculate the residues. Success with the method explained 
in Sect. 3.1 is predicated upon the assumption that the colour-singlet inhomoge- 
neous vertex can be described by a positive definite spectral density; e.g., this 
validates the form of Eq. (11). The method may fail, therefore, if the interaction 
produces a vertex that is inconsistent with this physical requirement. 

To explore that possibility we employed matrix inversion methods to solve 
Eq. (2) directly for timelike total momenta. The result, presented in Fig. 6, reveals 
a signal defect in the interaction model proposed in Ref. [37]. The model produces 
a vertex that possesses a zero at P 2 = —1.15 GeV 2 , and a zero in £^5(0, P 2 ) is 
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Figure 6. Solution of Eq. (2) obtained directly at timelike momenta. Three singularities 
that would normally indicate bound state poles are readily apparent. However, the circle at 
P 2 — —1.15 GeV 2 marks an unanticipated zero. 

only possible if the spectral density in the colour-singlet pseudoscalar channel 
is negative on a measurable domain. Similar behaviour is found in the scalar 
channel. 

An explanation of Figs. 4 and 5 is now evident. The zeros in 1/E^(0, P 2 ) 
provide a strong signal and may readily be recovered via a Pade approximant. 
However, it is plain that the singularity in 1/^(0, P 2 ), produced by the zero 
circled in Fig. 6 and lying just a little further into the timelike region than the 
position of the first excited state, can distort the value inferred for that state's 
residue. 

We have thus arrived via a circuitous route at important news. In attempting 
to obtain knowledge of bound states from information provided solely at spacelike 
momenta, we have uncovered a weakness in the interaction model of Ref. [37] that 
can plausibly be said to limit its domain of reliability to bound systems with 
mass < 1.1 GeV. However, in order to establish this we had to forgo the notion 
of using spacelike information alone. Without knowing the correct values of the 
residues a priori, one may reasonably have accepted the results in Figs. 4 and 5 as 
true measures of the first excited states' properties. The analysis in this section 
emphasises that in relying solely on spacelike information one rests heavily on 
assumptions made about the analytic properties of the Schwinger functions under 
examination. 

3.3 Configuration space analysis 



It is natural to ask whether more information about bound states can be ex- 
tracted from Schwinger functions specified in configuration space. These are, 
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for example, the quantities most often determined in numerical simulations of 
lattice-regularised QCD. A relevant comparison can be obtained with the results 
of Sect. 3.1 by analysing the Fourier transform of Eq. (11): 

M-l 

C(T) = W(T) + £^e-^\ (13) 

8=0 * 

where r is a Euclidean time variable. Equation (13) exhibits the behaviour ex- 
pected of a Schwinger function in a channel with only stable bound states. 

3.3.1 Effective mass 

For any observable projection of a two- or three-point Schwinger function the 
evolution at large Euclidean time is determined by the lowest mass state that 
couples to the channel under consideration. One can therefore define an effective 
mass whose large-r behaviour provides an approximation to the ground state 
mass: 

= m + o(e"( mi - m °) T ) . (14) 

The effective mass for the Schwinger function in Eq. (13) is depicted in 
Fig. 7. It was obtained from C(r) calculated at intervals St = 1 GeV -1 with 
r max = 16GeV~ 1 (> 1/mi). This is a mockup of the situation in which C(r) is 
only available on a discrete grid. The grid spacing used here does not affect the 
illustration in this subsection. As apparent in the figure, the rate at which m e s (r) 
approaches the ground state's mass depends on Am = mi — mo. 6 In principle, 
from a given Schwinger function the mass of the first excited state can be inferred 



m eff (r) :-- 



5t 



In 



C(t + St) 
C(r) 



6 The effective mass approaches a plateau from below because aj+i/tti < in Table 1. 
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Table 2. Masses (fJ,i, in GeV) extracted via a least-squares fit of a sum of N e exponentials, 
Eq. (15), to C(t) in Eq. (13) sampled at N = 16 intervals of St — 1 GeV - with r max = 
16GeV _1 . R 2 = 5^fL 1 [-EjVe( r /C{ri) — 1\ 2 . The last row repeats, for convenience of comparison, 
the generating values for C(r) (Table 1). 
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via the m e fj(r)-method if the ground state can be projected out. To be successful 
this requires that: the ground state residue and mass be accurately determined; 
and the level separation between the first and second excited states is not signif- 
icantly smaller than that between the ground and first excited state (a smaller 
separation requires that a larger domain of r have a good signal-to-noise ratio). 



3.3.2 Exponential fits 

The last observation completed logically suggests that one might recover the 
bound state spectrum contained in a given Schwinger function by fitting a sum 
of exponentials [43]: 

N e 

EnJj) =J> fc e-^ T (15) 
k=l 

Table 2 contains the results of a least-squares fitting procedure when applied to 
C(t) in Eq. (13), sampled as described above. 

This procedure is as efficient but not more so than that of Sect. 3.1. The 
lowest mass state is reliably obtained and, as N e is increased, so is the next 
least massive state. However, the remaining masses are false, even though the 
R 2 error is small. This fact is not changed by choosing instead an absolute- 
error measure for R 2 . This is disconcerting but straightforward to understand: 
for nonzero error tolerance a least-squares nonlinear minimisation problem does 
not generally possess a unique solution. Table 2 makes this plain empirically. 

The results in Table 2 can be viewed as an attempt to extract the masses 
within the framework of Bayesian statistics. The sample is fitted with a single 
exponential. The result is then used as an input constraint on a fit with two 
exponentials and so on, logically, until N e = 5. It is apparent that even with a 
subjective but accurate assumption about the masses of the two lightest states, 
nothing is gained in reliability of the estimates for the higher mass exponentials. 
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3.3.3 Error modelling 

The analysis described hitherto is based on a curve that is sampled discretely 
but with near perfect accuracy. This is the nature of DSE and kindred studies, 
in which the error is primarily systematic; i.e., it arises from truncation and, e.g., 
from uncertainties inherent in modelling the long-range part of the interaction 
between light-quarks. However, in principle, once truncation and model are cho- 
sen, machine precision is achievable. That is not so, e.g., with lattice simulations. 
In addition to systematic errors, such as those associated with finite size and 
volume, quenching, extrapolation to realistic current-quark masses, etc., there 
is a material statistical error. It is therefore interesting to generalise the model 
Schwinger function, Eq. (13), by adding Gaussian noise in the absolute value. 

To be specific, we generate M = 100 Gaussian-distributed random numbers 
{Q.j\ J = 1, . . . , M} with a mean of zero and a standard deviation of unity. Then 
from the exact values of C(r) sampled at N = 16 equally spaced points within 

[0, T max ] ■ 

{a = C{n) I n = ia,i = 1, . . . , N; a = r max /iV}, (16) 
we construct the sets Cj with elements 



(c< J ) i = c i + cg J 



1,...,N. 



(17) 



NB. Our reference value for the magnitude of the error is 

do = 0.01 , (18) 

which is ~ C(2Na)/C(0); i.e., the strength of the signal in the far infrared. 
We thus obtain M sets that are represented by the rows in the following array: 



(Ci% 
(C5)a, 



(CDn 



N 



(19) 



{Cm)n 



Each column in this array, labelled by i = 1,...,N, may also be considered 
as a set, each of which we will call a distinct configuration. The configuration 
characterised by a given value of i is an M-element set of measurements of C(jj) 
with Gaussian standard deviation ? around a mean value Ci. 

Reconsider now the effective mass defined in Eq. (14). From each row of 
Eq. (19), define 



m ; 



In 



1,...,N-1. 



(20) 



This provides M measurements of the effective mass for each time step; viz., the 
columns in 



■/?? 



i- 



ml, 



m 



m 



2- 



2< 



m 



N-i 
l 

N-l 



(21) 



rn 



Al- 



111 



Al- 



ii I 



N-l 
M 
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The average value of the effective mass at a given time step is 

1 M 

™ l = IlE^- ( 22 ) 
J=i 

The statistical error in fh 1 is determined via a jack-knife procedure, which 
yields a modification of the standard-deviation. It is designed to account for 
correlations in the statistical fluctuations between consecutive measurements of a 
given quantity. Here the quantity is fh 1 and we want to account for the possibility 
that if, e.g., fh 1 is above the true average then there is also a bias for fh %+1 to be 
above the average. If there are no such correlations in the errors then the method 
reproduces the standard deviation. 

The procedure is implemented as follows. For a given value of i one considers 
the sample {m'j, J = 1, . . . , M}. The element m\ is omitted and one defines 

1 M 

™y=w^i52 mlj,i=1, --- ,N ~ 1 - (23) 

J=2 

This is repeated, but with the J = 2 element omitted instead, so that m\ is 
thereby obtained. By stepping in this way through all M elements one arrives at 
a set {mV; J = 1, . . . , M} . The jack-knife estimate of the statistical error in fh 1 
is then 

J=i 

We arrive in this way at a model for a data set that could arise via sampling 
of a system whose evolution is described by an effective mass m e s{r); i.e., 

m cS = {fhi ±a^i, i = 1, . . . , N - 1} . (25) 

NB. While it is not immediately relevant, in the following we omit entries for 
which the error so calculated is greater than one-half of the mean value. 

Our approach to mocking up a Gaussian statistical error is well defined. It may 
nevertheless underestimate the full extent of measurement errors. For example, 
in a system with light-quarks one has for mesons [43] 

„ e -(»m-rr*)T (26) 

noise 

where m n is the pion mass, and this ratio is plainly small for meson masses 
not too much greater than m w . Moreover, one might find that in some channels 
there is a constant background noise, such as is the case in numerical simulations 
of glueball correlators [44]. This possible defect in our method will, however, 
only tend to make us underestimate the impact of errors on the reliability of an 
extraction of information from a data sample. 
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Figure 8. Effective mass plots for a system characterised by C(r) in Eq. (13) with Gaussian 
noise added (Sect. 3.3.3, M = 100 configurations): left panels, ground state; and right panels, 
first excited state, obtained with ao = in Eq. (13); i.e., with the lowest mass state omitted. 
Upper panels ~ c = 0.01 ?o, lower panels - ? = 0.1 <;o- In each panel the horizontal dotted line 
indicates the actual mass for the relevant state. 



3.3.4 Results from noisy data 

In Fig. 8 we illustrate the "noisy" effective mass generated from C(r) by the 
method just described. This figure may be compared with Fig. 7 and the impact 
of the error is readily apparent. 7 With q just 10% of <^o a mass plateau is quite 
effectively obscured. We subsequently work most often with q = 0.01?o- In this 
case we infer an effective mass by supposing that the domain [rg,ri5] lies within 
a mass plateau. This yields mass values that are too low by only 1%. 

We also reconsider the method of Sect. 3.3.2 in connection with the noisy 
sample. In this case one performs an N e exponential least-squares fit to each 
one of the M configurations in Eq. (19) and thereby determines the masses (and 
residues) for that configuration: {/i^, k = 0, . . . , (N e — 1)}. The mean fit-masses, 
/2fc, are evaluated via Eq. (22) and the jack-knife procedure employed thereafter 
to determine the statistical error in each of these. Results obtained for fio,fii via 
iV e = 3,4 exponential fits are depicted in Fig. 9. Within errors, they are consis- 

7 It is noteworthy that our error model produces data which has the appearance of that obtained 
in numerical simulations of lattice-QCD (see, e.g., Refs. [45, 46, 47, 48]). The magnitude of that 
error lies between our ? = 0.01 ?o and ? = 0.1 ?o samples. 
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Figure 9. r m i n -dependence in a least-squares fit of Eq. (15) to C(r) in Eq. (13) specified on 
t e [r m in,ri6] with Gaussian noise added (Sect. 3.3.3: M = 100 configurations; <; = 0.01 Co; 
R 2 < 10 -5 , uniformly). Left panel - N e = 3; right panel - N e — 4. Plateaux are visible in both 
cases. In both panels the horizontal dotted lines indicate the relevant mass values, m or mi. 



tent with those reported in Table 2. This figure contains additional information; 
namely, it illustrates the dependence of the extracted mass on the domain of r for 
which information is available. We performed an exponential fit on the domain 
t £ [T m i n ,T max = Tie], with T m i n = T\ , and then repeatedly incremented r m i n 
until the number of fit parameters exceeded the information sample. This occurs 
at r m i n = rii for N e = 3 (six pieces of information to fit six parameters) and 
Tmin = tq for N e = 4. The value of r m i n does not strongly influence the extracted 
value of ftQ. However, the error in p,\ tends to grow with T nnn . This is not surpris- 
ing because the ground state comes increasingly to dominate the signal at larger 
values of r. 

These two examples suggest that the analysis of data with a Gaussian error 
of <; = 0.01 ?o will yield results whose reliability is not significantly poorer than 
that obtained from a noiseless sample. Naturally, an improved result cannot be 
obtained from noisier data. Nevertheless, the techniques we have hitherto illus- 
trated are patently inadequate to the task of inferring reliably anything but the 
lightest two masses. A more advanced tool is necessary. 



4 Masses from Schwinger Functions: Correlator Matrix Method 

A more sophisticated method follows from an appreciation that the simultaneous 
analysis of Schwinger functions obtained from numerous carefully chosen interpo- 
lating fields can assist in determining the properties of states other than the one 
with the lowest mass in a given channel [32, 33]. This procedure will naturally 
become very time consuming as the number of operators is increased. 

To introduce the method [31, 47, 48], suppose one has operators {xi(x, r),i = 
1, . . . , K} that each have a nonzero residue in a particular channel and for which 
one has available the correlation matrix 

Gy(r,p) = J d 3 xe-^ x (01x^,^(0,0)10) . (27) 
Now insert a supposedly complete set of colour-singlet states whose quantum 
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£ \H n )(H n \ = l, 



(28) 



n=0 



where momentum and other labels are not written explicitly, and it may be that 
-^states = oo. One then has, using translational invariance of the vacuum, 



G ij (r):=G ij (r,p = 0)= £ e~ m ^ 



(29) 



n=0 



with 



^ = <0|x<(O,0)|fl n ), ^ = (iJ- n |xi(O,0)|0). (30) 

Equation (29) is a many-operator analogue of Eq. (13) in Sect. 3.3. NB. In the 
real- world application of this method, such as in numerical simulations of lattice- 
QCD, considerable effort must be devoted to the construction of the operators 
Xi in order to optimise for success. The operators are typically combinations of 
nonlocal few-quark/antiquark operators built and tuned empirically so as to have 
large overlap with states of interest in a given channel [31, 47, 48]. Our analysis 
will highlight the importance of this step. 

The correlation matrix G(r) is Hermitian and hence one can solve the eigen- 
value problem 



G{t q )- 1 ' 2 G{t)G{t ) 



'1/2 



<p k (r,r ) = \ k (r,T )4> k (r,ro) 



(31) 



to obtain K "principal correlators" A fc (r, to). 8 With the operators Xi defined and 
ordered such that X h > X k+l , it can then be shown that 



lim A fe (r,r ) = e - (r - To)m ( 1 + 0(e~ TAm ) 

One therefore has K principal effective masses 

-\ 1+1 {t + 5t,t q )- 



mm 



m J —m 



hi 



l = 0,...,(K-l) 



(32) 



(33) 



each of which will, for large r, exhibit a single plateau that corresponds to the 
mass of one of the K lowest-mass states in the channel under consideration. In 
practice K < Abates • 

We now explore and illustrate the efficacy of this method. Suppose there is a 
colour-singlet channel that possesses a trajectory of bound-states whose masses 
are given by [49]: 



mf = ml + Ifi 2 ,1 = 0,1, ... , (l T 



1), 



(34) 



where a T rriQ = 0.5, a 2 ./! 2 = 0.36 and the length-scale a T = 0.1 fm. We will subse- 
quently work with Z max = 20. 



Vo is a reference time that is typically assigned a small value. Reliable results should be inde- 
pendent of tq. We use ro = 0. 
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4.1 Perfect operators 

As a first step we suppose that a correlator matrix of the form in Eq. (29) has 
been obtained with K = 6 operators, and assume in addition that these operators 
are perfect; namely, each operator produces a unique state from the vacuum, so 
that 



The correlator matrix is therefore 



5 i(ra+l) 



(35) 



G(r) 



V 



— mor 




















e -miT 




















e -m 2 T 




















e -m 3 r 




















e -7Tl4T 




















e -m 5 r 



J 



(36) 



This situation cannot be achieved in reality because it requires a priori complete 
information about all bound states in a given channel and presumes that all the 
states are stable. It is nevertheless a useful demonstration case. 

To this correlator matrix we add Gaussian error in the manner described 
in Sect. 3.3.3; namely, the matrix is treated in the same way as the function in 
Eq. (13). A given element is sampled at N = 16 time steps, measured here in 
units of a T , to produce {Gij(sa T ), s = 1, ...,N}. We then generate M = 100 
Gaussian-distributed random numbers {Qj; J = 1, . . . , M} with a mean of zero 
and a standard deviation of unity, and thereafter construct the sets {Gij(s)Yj 
with elements 

[GijWj = Gij(s) + ? Gj ; a = 1, . . . , N . (37) 

This yields M correlator matrices sampled at N r-values, which appear as the 
rows in the following matrix: 



G(l){, G(2)\, 



cm 



G(i; 



G{2% 



G{2f 



G(N)\, 
G(N)l 



G(NY 



(38) 



(NB. Each entry is a K x K matrix.) As before, each column in this array, 
labelled by s = 1, . . . , N, may also be considered as a set, each of which is a 
distinct configuration. The configuration identified by a given value of s is an 
M-element set of measurements of the correlator matrix G(s) with Gaussian 
standard deviation <; around a mean value G(s). 

From each element in a given row of the matrix in Eq. (38), the K principal 
correlators can be determined by solving Eq. (31), and therefrom the principal 
effective masses: 



m l j(s) 



In 



Xf\8 + l,8o) 



\ l j +1 (s,s ) 



l = 0,...,(K -1) 



(39) 
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Figure 10. Effective mass plots from analysis of the correlator matrix in Eq. (36) with Gaussian 
noise added via M = 100 configurations with ? = 0.01 ?o. In each panel the horizontal dotted- 
line indicates the input mass, Eq. (34). 



This provides M measurements of each of the principal effective masses at each 
time step r(= sa T ); i.e., the columns in 



m[(l), mi(2), mi(s-l) 
m l 2 (l), m l 2 (2), m l 2 (s-l) 



l = 0,...,(K-l), (40) 



"4/(!)> m l M (2), m l M (s-l) 
and hence an average value of a principal effective mass at a given time step: 



^ {s) M 

j=i 



M 

rrE 7 ^^- ( 41 ) 



The statistical error in m l (s) is calculated via the jack-knife procedure described 
in Sect. 3.3.3 

The results of this analysis are presented in Fig. 10. With perfect operators 
a plateau at the correct mass [input, Eq. (34)] is immediately apparent in each 
panel, a feature that illustrates the power of the method. The signal is degraded 
through noise as r increases but the impact of that is minimal when the plateaux 
are so readily evident. 

It will be observed that in each panel some time steps are omitted and, for 
all but the ground state, effective mass values are not plotted for r greater than 
some minimum value. The omission occurs because Gaussian noise can lead to 
X(s + 1)/A(s) > 1 in which case Eq. (39) cannot define a positive mass. This fact 
and our decision to omit data for which the error is greater than one-half of the 
mean value explains the absence of points at larger r. 9 



9 Such features are also apparent in the data obtained from numerical simulations of lattice- 
QCD, e.g., Refs. [47, 48]. 
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4.2 Untuned cf. optimised operators 

It is impossible in principle to begin with perfect operators, Xi- Thus to simulate 
a more realistic situation involving the determination of a three-point Schwinger 
function we write 

e? = (-l) n \Ti\,i = l,...,K, (42) 

where "rj" is a number chosen randomly from a normal distribution with mean 
zero and variance one. 10 We subsequently normalise such that 

JVstates — 1 

\£f\ 2 = l,i = l,...,K, (43) 

n=0 

where in our concrete examples Abates = 'max- Normalisation does not compro- 
mise generality and is a useful precaution against precocious noise. Equations 
(42), (43) specify a rank K = 6 correlator matrix associated with a channel in 
some system that contains £ max = 20 bound states. To this matrix we added 
Gaussian noise, as in the case of perfect operators, and repeated the analysis 
that led to Fig. 10. 

The results are presented in Fig. 11. It is evident from column-one that the 
method is far less reliable with untuned operators since even a perfect measure- 
ment of the correlator matrix is unable to provide an accurate result for the third 
excited state: there is no true plateau. Furthermore, for this state the presence of 
even an extremely small amount of Gaussian noise in the correlator matrix elim- 
inates the possibility of obtaining any information: the signal is lost altogether 
with a Gaussian error of only <; = 0.1% so = 0.001 To = 0.00001. 

Our illustrations have shown that tuning the operators is an important step 
in the practical application of the correlator matrix method. This operation may 
be represented by emphasising the diagonal elements in the correlator matrix but 
allowing for leakage into the off-diagonal elements. To be specific, we write 



I? 



Vy, i = n + l 



i-v , , x1 ( 44 ) 



(-^states — 1) 



where again "r" is a number chosen randomly from a normal distribution with 
mean zero and variance one, and subsequently normalise. A value of y = 0.5 is 
described as 50% optimisation. This procedure specifies a correlator matrix to 
which we add Gaussian noise. 

The results for the principal effective masses are depicted in Fig. 12. With 
50% optimisation a clear plateau is apparent for the lowest mass state and, ar- 
guably, for the next lightest state, but beyond that one cannot confidently locate 
a plateau. Moreover, if one nevertheless chooses to estimate masses from defined 
rather than evident plateaux in the data, those masses are degenerate within er- 
rors. However, for 70% or better optimisation, clear plateaux are distinguishable 
for all of the six lightest states and accurate masses can be inferred. 



10 As noted in Ref. [17], the coefficients of states in a physical three-point Schwinger function 
must alternate in sign. This is not the case for a four-point polarisation function. 
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Figure 11. Principal effective masses from analysis of the untuned correlator matrix described 
in connection with Eqs. (42), (43). In this case the operators i = 1, ■ • ■ , K} are not perfect. 
Column 1 - Precise measurement of the correlator matrix; i.e., Gaussian noise ? = 0; Column 2 
- Gaussian error of c. = 1CP J qj; and Column 3 - Gaussian error of c = a. The horizontal 
dotted-line in the top nine panels indicates the relevant input mass; viz., mo, 1,2 from Eq. (34). 
In each of the bottom three panels there are two horizontal lines: the lower is 7713, which is the 
correct result, and the upper is 7714 in order to provide a context for any putative plateau. 



4.3 Influence of level location 

Two other features of the spectrum in a given channel may reasonably be expected 
to impact upon the extraction of masses from a Schwinger function; namely, the 
position of the ground state and the level separation. One may investigate these 
effects by varying the parameters in Eq. (34). To explore the first, in Fig. 13 we 
compare a T rriQ = 0.071 (mo ^ 0.14 GeV) with a T mo = 1.02 (mo — 2.0 GeV), while 
in Fig. 14 we contrast a^/i 2 = 0.18 with a 2 .// 2 = 1.02. NB. For these comparisons 
we used a tuned correlator matrix constructed according to Eq. (44) with 80% 
optimisation. 

It is evident in Fig. 13 that the primary effect of an increase in a r mo is a 
reduction, for all states, in the domain on which a plateau may be observed. 
The larger value of a r mo causes the signal to more rapidly enter the domain 
on which noise becomes dominant. Figure 14 illustrates that similar statements 
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Figure 12. Principal effective masses for the six lightest states generated by Eq. (34) calculated 
from the tuned-operator correlator matrix described in connection with Eq. (44) (Gaussian 
noise: M = 100 configurations; ? = 0.01 ?o): Upper 6 panels - 50% optimisation; Middle 6 - 
70%; Bottom 6 - 90%. The horizontal dotted-line indicates the relevant input mass. Where 
it appears, the solid line is an error-weighted least-squares fit to the data beginning with the 
point at r = 3a T . 



apply when changes in o? T \J? are considered. It follows that for any given channel, 
ground state mass and level spacing, smaller values of a T will help improve the 
method's accuracy if mnT max , r max = Na T , can be held fixed and noise levels 
kept constant. Satisfying these requirements simultaneously is challenging. 
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Figure 13. Principal effective masses obtained from a correlator matrix for the six lightest 
states generated by the formula in Eq. (34): upper 6 panels - a T mo = 0.071; lower 6 panels - 
a T mo = 1.02. In all panels the horizontal dotted-line indicates the relevant input mass. Where it 
appears, the solid line is an error-weighted least-squares fit to the data beginning with the point 
at r = 3a T . (Tuned correlator matrix with 80% overlap, Eq. (44). Gaussian noise: M = 100 
configurations; ? = 0.01 ?o). 



5 Epilogue 

We analysed the capacity of Schwinger functions to yield information about 
bound states, and established that for the ground state in a given channel 
the mass and residue are accessible via rudimentary methods. When the mass- 
splitting between the first and second excited states is not significantly smaller 
than that between the ground and first excited states, these methods may also 
reliably provide the same information about the first excited state, so long as 
the Schwinger function under consideration does actually possess the analytic 
properties assumed of it. However, simple methods cannot provide dependable 
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Figure 14. Principal effective masses obtained from a correlator matrix for the six lightest 
states generated by the formula in Eq. (34): upper 6 panels - a^/J. 2 = 0.18; lower 6 panels - 
a 2 -^ 2 — 1.02. In all panels the horizontal dotted-line indicates the relevant input mass. Where it 
appears, the solid line is an error-weighted least-squares fit to the data beginning with the point 
at r = 3a T . (Tuned correlator matrix with 80% overlap, Eq. (44). Gaussian noise: M = 100 
configurations; ? = 0.01 ?o). 



information about more massive states in a given channel. 

Indeed, there is no easy way to extract such information. An approach based 
on a correlator matrix can be successful but only if the operators are carefully 
constructed so as to have large overlap with states of interest in a given channel 
and statistical errors can be made small; viz., ~ 1%. While it is possible in princi- 
ple to satisfy these constraints, doing so is labour intensive and time consuming. 
That is nevertheless justified in the absence of model-dependence. 

We have posed the question of whether, in the context of bound-state studies 
in which model assumptions are made regarding the nature of the long-range 
interaction between light quarks, labour can be saved and/or accuracy gained by 
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working solely with Schwinger functions. In our view the analysis herein suggests 
strongly that the answer is no. In their formulation the studies of this type 
explicitly define an analytic continuation of a model into the complex plane. 
Hence all momenta are directly accessible. We saw clearly that it is not possible to 
avoid model artefacts by pretending ignorance of this. Hence nothing is gained by 
solving a complex of DSEs for the Schwinger functions associated with numerous 
carefully tuned interpolating fields. Moreover, capitalising on the feature that all 
momenta are directly accessible is the most efficient manner by which to test 
and improve the model input, and thereby to relate experimental data on bound 
state properties to information about the interaction between light quarks. This 
is further emphasised by the fact that because excited states are unstable they 
are characterised by complex pole positions. 
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Appendix A: Pole residue in a vertex 

Equation (2) is the inhomogeneous pseudoscalar Bethe-Salpeter equation (BSE). The general 
solution has the form in Eq. (4). Each scalar function therein depends on the renormalisation 
point but we have not written that explicitly. The inhomogeneous equation is solved subject to 
a renormalisation boundary condition: 

which translates into the single condition 

E 5 (k-P)\ k2=(2p=0 = l. (A.2) 

This result is guaranteed so long as the renormalisation constant for the pseudoscalar vertex is 
the same as that for the Lagrangian mass; viz., 

Z 5 (( 2 ,A 2 ) = Z4{ 2 ,A 2 ). (A.3) 

This is a consequence of the axial-vector Ward-Takahashi identity. 
In the neighbourhood of any pseudoscalar meson pole, 

r i^ P )\ 2 , = P 2+" 2 r L(k;P), (A.4) 

where 

ri n (k; P) = r J 75 [iE nn (k; P)+y PF Vn (k; P) + 7 • k k ■ P G* n (k; P) + k»P v H Wn (k; P)] 

(A.5) 

is the canonically normalised bound-state Bethe-Salpeter amplitude, and 

ip^iC 2 ) = MC 2 , A 2 ) tr f 75 Xx n (q; P) (A.6) 

is the gauge-invariant and regularisation-scale-independent pseudoscalar projection of the me- 
son's Bethe-Salpeter wave function onto the origin in configuration space. In Eq. (A.6), x^n is 
specified analogously to Eq. (3) but with the inhomogeneous vertex replaced by Eq. (A.5). Here, 
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as elsewhere herein, the suffix n labels nondegenerate states in the channel: n = is the ground 
state, n — 1 is the first excited state, etc. 

In the isospin symmetric limit, Eq. (A. 4) entails 

tr F r'rUk; P) = S ij r 5 (k; P) = 5 ij ^ " - P^ n (fc; P) (A.7) 

so that in the neighbourhood of a pole 

r 5 (0;P)| p2+m ,^ = p2 2 ^^ n (0;P) (A.8) 

or, equivalently, 

g s(0; P)= nTT 2 2p„ n E„ n {0;P). (A.9) 

The residue of the pole in the vertex is therefore 

r E „ n = 2p„„ E„„ (0; P) , (A.10) 

which can be calculated so long as one has the solution of the homogeneous BSE in addition 
to the solution of Eq. (2). Thus do we calculate the dotted curves in the right panels of Figs. 4 
and 5. 
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